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Abstract. Anuran amphibians have intensively been studied to understand Amazonian biodiversity. Improved methods 
and sampling has revealed that many widespread nominal species in fact are complexes of species with smaller allopatric 
ranges. Pan-Amazonian anuran species are rather an exception. In a case study using the three-striped poison frog (Anura: 
Dendrobatidae: Ameerega trivittata), we ask how the pan-Amazonian distribution of this taxon can be explained and hy- 
pothesize that dispersal has played a major role. Species delimitation and intraspecific relationships of the study species 
were examined from novel and existing (GenBank) sequences of the mitochondrial 16S rRNA gene from 108 specimens of 
38 localities using maximum likelihood and Bayesian methods. We performed BioGeoBEARS models using a time-cali- 
brated population tree to reconstruct the biogeographic history. Our results support that A. trivittata is a pan-Amazonian 
species scattered over its geographic range. Being of Late Miocene origin, the species rapidly spread into newly available 
space and repeatedly dispersed for-and backward, while vicariance played a major role only in the Early Pliocene. We sug- 
gest that intrinsic morphological and life history characteristics (adult size, relative reproductive success) make A. trivitta- 
ta a more successful disperser than other species, so that riverine barriers are more permeable and hamper allopatric 
speciation. We conclude that there is no universal causality explaining Amazonia biodiversity, because species-specific 
biological characteristics are key determents of biogeographical histories. Comparatively better dispersal advantages foster 
larger geographic ranges and can explain pan-Amazonian distributions. 


Key words. Amphibia, Anura, BioGeoBEARS, dispersal, genetic diversity, mitochondrial DNA, Neotropics, riverine bar- 
rier, vicariance, widespread species. 


Introduction 


The Amazon basin of South America is one of the regions 
with the highest biological diversity on earth (JENKINS et 
al. 2013, ANTONELLI et al. 2015). There have been various 
hypotheses about the diversification processes of its extant 
biota (HILL & HILL 2001, Hoorn et al. 2010, RULL & CAR- 
NAVAL 2020). However, studies on both plants and animals 
revealed that there is no universal causality to explain spe- 
cies richness in Amazonia. Rather, for particular groups, 
the biogeographical history is more unique. Reasons in- 
clude species-specific life history traits, biotic interaction, 


habitat availability, river course change, geological age and 
geographic origin (COLLEVATTI et al. 2009, LEITE & RoG- 
ERS 2013, ANTONELLI et al. 2018, CRACRAFT et al. 2020,). 
This is well exemplified in the more than 600 extant am- 
phibian species, most of which belonging to the order An- 
ura (e.g. MAYER et al. 2019). A variety of hypotheses have 
been invoked to explain the species richness within differ- 
ent taxa (reviewed by LEITE & ROGERS 2013). 

Some examples for anuran taxa are given in the fol- 
lowing. NOONAN & GAUCHER (2005) demonstrated in- 
complete speciation processes due to repeated secondary 
contact in harlequin toads of the genus Atelopus from the 
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eastern Guiana Shield in the frame of the Disturbance Vi- 
cariance hypothesis. The widely invoked river-barrier hy- 
pothesis was demonstrated to be applicable to tree frogs 
(Dendropsophus leucophyllatus complex) by Prrani et al. 
(2019), suggesting dispersal rather than vicariance scenar- 
ios. REJAUD et al. (2020), using poison frogs of the genus 
Allobates, showed that allopatric speciation (i.e. vicariance) 
played a key role for explaining species richness in the up- 
per Amazon basin, while the lower basin and its vicinities 
were colonized mainly via dispersal. 

So far, there is no hypothesis explaining the evolution- 
ary biogeography of pan-Amazonian anurans. There might 
be a simple reason for this. While in the past, many Ama- 
zonian anurans were thought to encompass large distribu- 
tions - ranging from the Andean versant in the west upper 
basin to the Guiana Shield in the north-east (e.g. LYNCH 
1979) — more recent studies demonstrated that this is not 
the case. Many of the nominal species with large Amazoni- 
an distributions actually represent complexes of allopatric 
cryptic taxa with much smaller geographic ranges (Fou- 
QUET et al. 2007, 2012, 2014, PADIAL & DE LA RIVA 2009, 
ANGULO & ICOCHEA 2010, BROWN et al. 2011, FUNK et al. 


2012, GEHARA et al. 2014, PELOsO et al. 2014, VACHER et 
al. 2017, Rojas et al. 2018). As a result, truly pan-Amazo- 
nian anuran species are rather an exception that yet have 
received little attention. 

One pan-Amazonian species is the three-striped poi- 
son frog Ameerega trivittata (SP1x, 1824), a member of the 
family Dendrobatidae, subfamily Colostethinae (Fig. 1). It 
is one of the largest (adult snout-vent length 35-55 mm) 
and most widespread species of the entire superfamily 
Dendrobatoidea (GRANT et al. 2006, LOTTERs et al. 2007). 
The conspecificy of morphs allocable to A. trivittata from 
different Amazonian terra firme forest localities has been 
proposed on the basis of morphology (SILVERSTONE 1976) 
and molecular genetics (ROBERTS et al. 2006, GRANT et 
al. 2017, GUILLORY et al. 2020). According to SANTos et al. 
(2009), A. trivittata originated in the Amazon basin in the 
Late Miocene and GUILLORY et al. (2020) suggested that 
this was around 7.57 Mya. 

Based on literature information on the geological his- 
tory of the Amazon basin, we develop some scenarios 
that may help explaining the evolutionary biogeography 
of A. trivittata. With the emergence of the Amazon run- 


Figure 1. Schematic map of northern South America showing localities of Ameerega trivittata investigated in this study (white sym- 
bols) and their allocation to biogeographic regions (BR1-BR7 after GODINHO & SILVA 2018, including colors). Samples of the Madeira 
River group (see text) are shown as squares, whereas all other samples are given as dots. Two localities outside BRs were provisionally 
allocated to their nearest BR (indicated by red arrows; Soldado Oliva and Puerto Maldonado, respectively, Table 1). The geographic 
range of A. trivittata as suggested by the IUCN is indicated by a solid line (https://www.iucnredlist.org, accessed 13 February 2022). 
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ning into the Atlantic Ocean and starting around 8 Mya 
(SHEPHARD et al. 2010), A. trivittata may have expanded 
its distribution within central Amazonia. Such a scenario 
is postulated in various taxa (Hoorn et al. 2010) and we 
assume this for A. trivittata, too. Subsequent to that, start- 
ing circa 7-6 Mya, enormous water retention bodies for- 
matted in western and central Amazonia, that today are 
known as major Amazon tributaries (cf. HOORN et al. 2010, 
SHEPHARD et al. 2010). It is highly expectable that these 
water bodies were barriers to A. trivittata (as a terra firme 
species), so that populations became separated. However, 
when around 5 Mya water bodies opened into the Ama- 
zon and became rivers (HOORN et al. 2010, ALBERT et al. 
2018), they were perhaps permeable to A. trivittata allow- 
ing multi-directional dispersal and genetic exchange. We 
base this assumption on the observation that, in anurans, 
larger body size has been shown to be beneficial for disper- 
sal (WOLLENBERG et al. 2011, SEARCY et al. 2018). 

Using a molecular genetic approach (1) we here test 
for conspecificy of pan-Amazonian A. trivittata popula- 
tions using an enlarged geographic sampling. We expect 
all samples to belong to a single species based on (A) low 
genetic diversity and (B) representing a monophyletic as- 
semblage. In line with the one-species-hypothesis, we ex- 
pect (C) that A. trivittata recently invaded into various 
directions. Examining the phylogenetic architecture in a 
spatio-temporal context and using biogeographic model- 
ling, (2) we test for and expect early (Late Miocene) vicar- 
iance and subsequent (Pliocene to Pleistocene) dispersal 
events in this species. 


Methods 
Sampling and data preparation 


Ameerega trivittata samples were collected from 2014 
to 2016 at different localities in Brazil and Peru (Tab. 1). 
For these collections with use the following permits: SIS- 
BIO45202-1/2, 15PE 000174/SP, (authorized by permits 
0196-2014-MINAGRI-DGFFS/DGEFFS, 0050-2015-SER- 
FOR-DDGGSPFFS-DGSPFS, 312-2015-SERFOR-DGGSPF 
FS-DGSPFS) and completed with samples from scientific 
collections. Samples were sequenced (see below). Addi- 
tional sequences were obtained from GenBank (BENSON et 
al. 2015; https://www.ncbi.nlm.nih.gov/genbank/, accessed 
18 December 2021). In total, the complete dataset com- 
prised 108 specimens from 38 localities throughout most 
of the known geographic range of A. trivittata (cf. SILVER- 
STONE 1976, LOTTERs et al. 2007). See Supplementary Ta- 
ble Si for sampling locality details and GenBank accession 
numbers. 

Of these, 36 localities fall into six of the seven biogeo- 
graphic regions (BR) suggested for Amazonian anuran di- 
versity (GODINHO & SILVA 2018). We provisionally allocat- 
ed the two remaining localities to the existing BRs in their 
vicinities (Fig. 1). In the manner of similar studies (e.g. RÉ- 
JAUD et al. 2020), we use these widely acknowledged spatial 
units for biogeographic modelling (see below). 


We conducted direct sequencing with PCR primers 
on an automated Sanger sequencer to obtain DNA se- 
quences of the mitochondrial 16S rRNA gene fragment. 
This marker was chosen because it is suggested to ade- 
quately reflect genetic variability within and between spe- 
cies (e.g. VENCES et al. 2007). In addition, this marker 
has recently been successfully used to study the biogeo- 
graphy of several Amazonian amphibians (e.g. RÉJAUD et 
al. 2020, FouQUET et al. 2021a,b). We analyzed 83 sam- 
ples of A. trivittata. Genomic DNA was isolated using the 
QIAGEN DNeasy Tissue Kit, following the manufactur- 
er’s guidelines. We performed polymerase chain reaction 
(PCR) with primers of the mitochondrial gene fragment 
16S rRNA (16Sal 5’-CGC CTG TTT ATC AAA AAC AT-3’ 
and 16Sbh 5’-CCG GTC TGA ACT CAG ATC ACG T-3’) 
of PALUMBI et al. (2002). This fragment was amplified us- 
ing an initial denaturation step of 2 min at 94°C, followed 
by 35 cycles of denaturation at 94°C for 30 s, annealing at 
54°C for 90 s and extension at 65°C for 90 s. PCR prod- 
ucts were purified using spin columns (QIAGEN). Stand- 
ard Sanger sequencing was performed directly using 
the corresponding PCR primers and ABI Big Dye on an 
ABI3730XL sequencer. Each sample was sequenced twice 
in forward and reverse directions and complementary se- 
quences were aligned (in total 611 bp). Chromatograms 
were checked and edited visually. Sequences were depos- 
ited in GenBank (for accession numbers see Supplemen- 
tary Table S1). 

We used Mega 10.1 (KUMAR et al. 2018) for sequence 
alignment and analysis. We trimmed initial and end por- 
tions of the sequences to 484 bp, so no missing data ap- 
peared. Sequences were aligned with Muscte with de- 
fault parameters (EDGAR 2004) as implemented in Mega. 
Haplotypes were identified and analyzed with DnaSP 
6.12.03 (Rozas et al. 2017). We created the final haplo- 
type alignment with GBlocks 0.91.1 (CASTRESANA 2000, 
TALAVERA & CASTRESANA 2007; https://ngphylogeny.fr, 
accessed 24 April 2022) using standard settings and al- 
lowed gap positions ‘with half’ (Massana et al. 2004, 
MANICHANH et al. 2008). This collapsed our dataset into 
36 haplotypes, in part occurring at multiple sample sites, 
resulting in 57 sequences which were used for all analyses 
(Tab. 1). 

TPMa2uf+I+G was identified as the best-fitting sub- 
stitution model with JModelTest 2.1.10 (DARRIBA et al. 
2012), based on the corrected Akaike Information Crite- 
rion (AICc). When this model was not implemented in the 
software used, we applied the K80+I+G or the GTR+I+G 
model (see below). 


Species delimitation 


Goal (1A) was to test if genetic samples assigned to A. tri- 
vittata represent the same taxon. For this purpose, we 
used three approaches: Pairwise distances (i.e. uncorrect- 
ed p-values), Generalized Mixed Yule Coalescent mod- 
el (GMYC; Pons et al. 2006, FUJISAWA & BARRACLOUGH 
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Table 1. Sampled Ameerega trivittata localities classified to biogeographic regions (BR) with identified haplotypes (H1-H36). Localities: 
* Madeira River group (see text); ** not in BRs, here provisionally allocated to BRs (see Fig. 1). 


Biogeographic Country Locality Longitude Latitude Haplotypes (N) 
region 
BRI Brazil Tefé -3.46444 -64.67972 H15(3) 
Brazil Alenquer -1.45639 -55.09750 H3(5) 
Brazil Lago Miriti -3.32161 -59.54158  H2(1), H20(1), H21(1), H26(1) 
Brazil Manaquiri -3.68948 -60.33507 H24(3), H25(1) 
Brazil Manaus 1 -3.35338 -59.85563 H22(1) 
Brazil Manaus 2 -3.66866 -60.30093 H2(2), H23(1) 
Brazil Castanho, 40 km S Manaus -3.61956 -60.45511 H2(3) 
Brazil Uatuma River, Balbina -1.91015 -59.54855 H3(1) 
BR2 Suriname Kabo Forestry Concession 5.26642 -55.75567 H3(1) 
Suriname Paramaribo-Apura road 5.55848 -55.42227 H3(1) 
Suriname Brokopondo 5.02026 -55.07899 H3(1) 
BR4 Peru Tahuayo -4.31400 -73.23259 H14(1) 
Peru Manati River -3.65201 -72.20045 H13(1) 
Peru Bolognesi -10.06923 -73.91740 _H27(6), H28(1), H29(1), H30(1) 
Peru Contamana -7.22536 -74.95925 H7(2), H12(5), H31(1) 
Peru Panguana -9.58333 -74.80000  H4(2), H32(3), H33(2), H34(1), H35(1), 
H36(1) 
Peru Sepahua -11.13100 -73.05800 H7(8), H27(1) 
Peru Pebas -3.31285 -71.85243 H6(1) 
Peru Shucushuyaco -6.03204 -75.85705 H9(1) 
Peru Cainarachi Valley -6.33289 -76.28336 H9(2) 
Peru Shapaja -6.58088 -76.26275 H9(1) 
Peru Chumilla -6.61859 -76.18432 H11(1) 
Peru Tarapoto -6.46161 -76.35080 H9(1), H10(1) 
Peru Cordillera Azul -6.93666 -76.23903 H12(1) 
Peru Soldado Oliva** -5.30416 -78.38802 H9(1) 
Colombia Leticia -4.12330 -69.94910 H6(3) 
BR5 Brazil Purus River, Camicua -8.72703 -67.42433 H5(4) 
Brazil Juruá River, Ipixuna -7.05225 -71.69469 H17(1) 
Brazil Juruá River, Seringal do Condor -6.74825 -70.78964 H18(1) 
Brazil Porto Walter -8.24875 -72.78001 H7(8), H8(1) 
Peru Puerto Maldonado** -12.83715 -69.29372 H5(2) 
BR6 Brazil Madeira River, Ilha do Búfalo* -9.14375 -64.51606 H1(3) 
Brazil Madeira River, Jirau Direito -9.32417 -64.71083 H5(1) 
Brazil Madeira River, Ilha da Pedra* -9.15889 -64.63389 H1(2), H16(1) 
Brazil Madeira River, Teotônio* -8.84788 -64.06847 H1(1) 
BR7 Brazil Madeira River, Castanho* -5.26545 -61.94445 H1(1) 
Brazil Madeira River, Humaitá* -6.56312 -62.93648 H19(2) 
Brazil Tapajós River, left side -5.06892 -56.87200 H21(2) 


2013), and Assemble Species by Automatic Partitioning 
(ASAP; PUILLANDRE et al. 2020). We calculated pairwise 
distances in Mega under the K80 model (100 bootstraps, 
pairwise deletion of gaps; KIMURA 1980). For GMYC, we 
built an ultrametric tree using BEAST 2.6.3 and BEAUTI 
for parameter settings (GTR+G+I as substitution model, 
a strict clock model, Yule model as prior). Analysis was 
run with 10,000,000 chains, whereof every 10,000" tree 
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was stored. We set the burn-in to 25,000. The resulting 
10,001 trees were further compiled with TreeAnnotator 
2.6.3 (BOUCKAERT et al. 2019) to a maximum clade cred- 
ibility tree with common ancestor heights by defining a 
burn-in of 10% and a posterior probability limit of 0.9. This 
final tree was then used for GMYC (single threshold ver- 
sion) using the R packages Paran (DINNO 2009; software 
retrieved from http://cran.r-project.org/web/, accessed 
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30 July 2022), Splits (EZARD et al. 2009; http://R-Forge.R- 
project.org/projects/splits/, accessed 30 July 2022), Mass 
(RIpPLey et al. 2015; software retrieved from http://cran.r- 
project.org/web/, accessed 30 July 2022), and Ape 5.0 
(PARADIS & SCHLIEP 2019; software retrieved from http:// 
cran.r-project.org/web/, accessed 30 July 2022). As a third 
species delimitation method, we created species partitions 
from single locus sequence alignments on the ASAP web 
server (https://bioinfo.mnhn.fr/abi/public/asap/asapweb. 
html, accessed 24 June 2021), with 'P' between 0.001 and 
0.1; X = 0.5; N = 10 and K8o as substitution model, with 
ts/tv = 11.55 (transition/transversion ratio, according to the 
AICc results; see below). 


Phylogenetic analyses 


To test for monophyly of A. trivittata samples (goal 1B), we 
performed Bayesian Inference (BI) and Maximum Like- 
lihood (ML) phylogenetic reconstructions. As outgroups 
we used Ameerega braccata, A. picta and Colostethus pratti, 
based on the topologies of GuiLLory et al. (2020). Out- 
group sequences were obtained from GenBank (acces- 
sion numbers: DQ502125, KJ940455, KR863143). We ran 
BI analyses in MrBayes 3.2.7 (RONQUIST & HUELSENBECK 
2003) with 100 million generations and two independ- 
ent runs starting from different random trees, using the 
GTR+I+G substitution model. Chains were sampled eve- 
ry 10,000 generations. Convergence between the two runs 
was checked in Tracer 1.7.1 (RAMBAUT et al. 2018). The first 
25% of the trees obtained were rejected as burn-in and the 
remaining sampled trees were used for building a consen- 
sus tree and estimation of Bayesian posterior probabili- 
ties (BPP). ML analyses were performed with IQtree 2.1.2 
(NGUYEN et al. 2015) under the TPM2uf+I+G model. Tree 
node support was assessed with 10,000 ultrafast bootstrap 
(UFBoot) replicates (MINH et al. 2013). For both analyses, 
base frequencies, proportion of invariable sites and gamma 
distributed rates were adopted according to the results of 
AICc. 


Haplotype diversity 


Goal (1C) was to test for signal of recent invasions into dif- 
ferent directions by studying haplotype diversity. Genetic 
diversity parameters were computed with DnaSP includ- 
ing the number of polymorphic sites (S), nucleotide di- 
versity (n), haplotype diversity (Hd) and mean number of 
nucleotide difference (K). Furthermore, we analyzed with 
DnaSP haplotype difference and the ratio of transitions 
and transversions. For intraspecific gene genealogies, we 
analyzed the phylogenetic relationships between sequenc- 
es by creation of a haplotype network as follows. A statisti- 
cal parsimony network was created in TCS 1.23 (POSADA & 
CRANDALL 2000), with the 98% criterion for a parsimoni- 
ous connection, then visualized with tcsBU (MUrIAs Dos 
SANTOS et al. 2016). 


Evolutionary biogeographical reconstruction 


To test for early vicariance events in the Late Miocene and 
subsequent dispersal of A. trivittata (goal 2), we used a dat- 
ed phylogeny and biogeographic modelling, similar to RÉ- 
JAUD et al. (2020). 

Time-Calibrated Population Phylogeny. - To estimate 
the intraspecific diversification through time, we conduct- 
ed a dated phylogeny using BEAST and BEAUTI for de- 
fining parameters (BOUCKAERT et al. 2019). There are no 
fossil records for dendrobatoids and their ancestors avail- 
able, which prevents the possibility of a primary calibra- 
tion of the molecular clock. We therefore applied a second- 
ary calibration derived from GuILLory et al. (2020). These 
authors estimated the split of Ameerega and Colostethus to 
21.797 Mya with o = 3.071 Mya, based on divergence time 
estimations with a fossil record background for the entire 
class of Amphibia as described in Santos et al. (2009). We 
used the estimated divergence times of the three-striped 
poison frog from its three close phylogenetic relatives 
(according to GUILLORY et al. 2020), Ameerega picta at 
7.567 Mya with o = 1.634 Mya, A. braccata at 9.252 Mya with 
o = 1.846 Mya, and Colostethus pratti at 23.325 Mya with o = 
3.071 Mya as calibration points. Again the TPM2uf+I+G 
substitution model was employed. Further substitution 
parameters were set to four gamma rate categories with 
a fixed value for gamma shape and a proportion invari- 
ant as defined by JModeltest. We set a relaxed log-normal 
clock model with a clock rate prior of 1e-10 and a Yule tree 
model with default settings, as endorsed by GUILLORY et al. 
(2020). We ran three independent chains of 100,000,000 
iterations with different random seeds. Trees were logged 
every 10,000 steps and the first 1,000,000 chains were ig- 
nored as pre-burn-in. ESS (effective sample size) values 
were controlled in Tracer and were above the recommend- 
ed threshold of 200. The three runs were combined with 
LogCombiner 2.6.3 (BOUCKAERT et al. 2019) and summa- 
rized with TreeAnnotater using a maximum clade credibil- 
ity tree and mean node heights with 10% burn-in. 

Modelling. - A reliable method that combines phyloge- 
netic data and geographic data to describe biogeographical 
events is the R package BioGeoBEARS of MATZKE (2013). 
We accomplished ancestral area reconstruction and bio- 
geographical event estimations under BI and ML frame- 
works in BioGeoBEARS. Here, the molecular clock tree 
obtained from BEAST served as an input. Our sampling 
contained various haplotypes from the same population, 
which can bias the results in BioGeoBEARS (N. MATZKE, 
in litt. 19 August 2022). We therefore carefully pruned our 
tree to a meaningful population tree, using the pruning 
skript “How (and whether) to collapse tips to prune a tree” 
on PhyloWiki (http://phylo.wikidot.com/example-biogeo- 
bears-scripts#pruning_a_tree, accessed 28 August 2022). 
Moreover, we followed the recommendations of recent 
studies (e.g. PIRANI et al. 2020, REJauD et al. 2020, Fou- 
QUET et al. 2021a, b) and compared AICc values of a disper- 
sal extinction cladogenesis model, representing vicariance 
and in situ diversification (DEC; cf. REE & SMITH 2008), a 
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likelihood version of the dispersal vicariance model which 
stands for a majority of vicariance events (DIVALIKE; cf. 
RONQUIST 1997), and a likelihood version of the BayArea 
model, preferred for the majority of in situ diversification 
(BBM; cf. Lanns et al. 2013). We also followed the cau- 
tionary advice of REE & SANMARTIN (2018) and did not 
include jump dispersal parameters, as these can enlarge 
cladogenetic events to the detriment of anagenetic events 
and time-dependent range evolution. Following Dupin et 
al. (2017), we conducted 100 independent runs of biogeo- 
graphical stochastic mapping (BSM) in BioGeoBEARS to 
determine biogeographical event counts for the best-fitting 
model. 

To finally reconstruct processes at different spatio-tem- 
poral scales, we allocated sample localities of A. trivitta- 
ta into the BRs (Fig. 1). According to GODINHO & SILVA 
(2018), these BRs well reflect the regionalization of anuran 
dissimilarity in Amazonia. The spatial clustering into BRs 
allowed us to investigate potential dispersal events of our 
study species during the diversification of the Amazon 
basin and the establishment of the Amazon River system 
from the Acre system over the last ca. 9 Ma (cf. VACHER 
et al. 2017). We allowed for the three models mentioned 
above a maximum of six reconstructed areas per node, 
which equals an occupation of all BR groups where A. tri- 
vittata occurs. 


Results 
Species delimitation, phylogeny and genetic diversity 


The comparison of the uncorrected p-distances and GMYC 
suggested within-species variation. In contrast, ASAP de- 
limited two groups suggesting the possible presence of two 
taxa under the name A. trivittata. In detail, the uncorrected 
p-distances revealed genetic variation at < 2.3%, with most 
values < 2% (Supplementary Table S2). Taking the rule- 
of-thumb threshold of 5% for anuran species delimitation 
in this fragment of the 16S rRNA gene (fide VENCEs et al. 
2005b), our A. trivittata samples rather represent conspe- 
cifics. Under the GMYC model, all A. trivittata haplotypes 
were delimited into three clusters (Supplementary Table 
S3). However, there was no geographic signal, and the like- 
lihood ratio test was not significant (LR test = 0.097), indi- 
cating that our A. trivittata dataset belongs to a single spe- 
cies. In contrast, the ASAP analysis revealed two groups, 
ie. potential taxa, with an ‘ASAP score’ of 4.00, a P-val 
rank of 0.93 and a threshold distance of 0.00729. These two 
groups comprised: (i) three haplotypes from five localities 
along the Madeira River: Castanho (H1), Humaita (H19), 
Ilha da Pedra (Hi, H16), Ilha do Bufalo (H1) and Teotônio 
(H1) (Figs 1, 2); (ii) all other samples including those from 
the type locality (Rio Tefé, Estado Amazonas, Brazil). 
Both phylogenetic trees (BI: Fig. Sı; ML: Fig. S2) sup- 
port the monophyly of A. trivittata. In concert with the 
ASAP results and the haplotype network, samples of the 
Madeira River group are defined as a sister clade to all oth- 
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Table 2. Genetic diversity of Ameerega trivittata haplotypes found 
in the respective biogeographic regions (BR) with N = number of 
samples; h = number of haplotypes; S = number of polymorphic 
sites; m = nucleotide diversity; Hd = haplotype diversity and K = 
mean number of differences. 


BR N h S II Hd K 

BRI 13 10 12 0.00848 0.949 4.103 
BR2 3 1 0 0 0.000 0.000 
BR4 27 14 13 0.00577 0.934 2.758 
BRS 6 4 4 0.00387 0.867 1.867 
BR6 5 3 0.00579 0.700 2.800 
BR7 3 3 8 0.01102 1.000 5.333 


all BRs 57 36 26 0.00829 


er A. trivittata with robust support (BI: BPP = 100; ML: 
UFBoot = 91). 

The 36 unique haplotypes (108 samples; Tab. 1) differed 
at 26 polymorphic sites with 10 singletons and 16 parsimo- 
ny informative sites and showed ts/tv = 11.55. Haplotype di- 
versity and nucleotide diversity of all identified haplotypes 
were relatively high (Hd = 0.966; Pi = 0.829%). Haplotype 
diversity varied across sample sites and biogeographic re- 
gions (Tab. 2). The maximum number of haplotypes at one 
locality was six (Panguana, at the same time the locality 
with the highest samples size); the most common haplo- 
type (H3) was recorded at five localities (Tab. 1). The haplo- 
type network shown in Figure 2a revealed a ‘star burst’ pat- 
tern for samples from western and central Amazonia (BR1, 
BR4, BRs). Genetically most distant from the star burst 
were the haplotypes from regions BR2, BR6 and BR7. In 
addition, the haplotypes of the Madeira River group (Hı, 
H16, Hig) from regions BR6 and BR7 formed a cluster that 
was not connected to the other haplotypes (Fig. 2a, b). 


Within-species evolution and biogeography 


In the independent BEAST runs, the A. trivittata samples 
were discriminated further (Fig. S3) resulting in (1) the 
Madeira River group, (2) a Northeastern central Amazo- 
nian group, (3) a Western central and southern Amazoni- 
an group. These splits were supported by BPP = 0.972 and 
0.902, respectively. Our time estimation dated the diversi- 
fication of A. trivittata from its sister A. picta at 8.44 Mya 
(95% HPD: 10.6-6.04 Mya). Ameerega trivittata split fur- 
ther at 6.87 Mya (95% HPD: 9.2-4.55 Mya), representing 
the basal split of the Madeira River group from groups (2) 
and (3). At 5.75 Mya (95% HPD: 7.91-3.56 Mya), the last two 
mentioned groups split. 

Comparisons of the BioGeoBEARS models AICc iden- 
tified DIVALIKE as the best-fit model (Tab. 3). Our area 
reconstruction revealed that the most likely ancestral es- 
timated area of A. trivittata was composed of BR1, BR4, 
BRs and BR6 (Fig. 3). Moreover, biogeographical stochas- 
tic mapping suggested that the segregation of the principal 


Biogeography of a pan-Amazonian poison frog 


groups (1) to (3) is best explained by vicariance events (V1 
around 6.87 Mya and V2 around 5.75 Mya; Fig. 3). BioGeo- 
BEARS suggested an additional vicariance event within 
group (3), which segregates a geographic unit comprised 
by (A) the southern Amazonian samples from (B) the west- 


ern central Amazonian samples (V3 around 5 Mya; Fig. 3). 
However, this split received no support in the BEAST to- 
pology (Fig. S3). Subsequent within-group diversification, 
younger than circa 5 Mya, was identified to be best ex- 
plained by dispersal events. 


Figure 2. Haplotypes of Ameerega trivittata: (a) Statistical parsimony network pattern (98%) inferred from the 36 haplotypes from 
38 localities, with the Madeira River group (H1, H16, H19, see text) showing a separate network (using a 95% parsimony limit: H1 
connects to H3, H5, H7 with five mutations; H19 connects to H6 with five; H16 connects to H6 via H19 with six mutations - not 
shown). Colors are according to BR allocation (Fig. 1, Tab. 1); the number of samples per haplotype is indicated by circle size on 
the left margin. (b) Map of northern South America showing the distribution of the unique and the most shared haplotypes (faded 


background colors are those of BRs in Fig. 1). 
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Table 3. BioGeoBears AICc results. Abbreviations: LnL = In(likelihood); d = dispersal; e = extinction; j = founder-event speciation. 


Model LnL d e j AICc AlCc wt 
DEC -63.7359 0.0152 1.0e-12 0 131.7718 5.2994 
DIVALIKE -59.0920 0.0207 1.0e-12 0 122.4839 5.5088 
BAYAREALIKE -80.5670 0.0112 0.2852 0 165.4340 2.5976 
Discussion an anurans that have been studied in detail in recent years 


Hypotheses testing 


Our results demonstrate low genetic divergence among 
(uncorrected p-distances, GYMC, but not ASAP, as dis- 
cussed below) and monophyly (BI and ML) of A. trivittata 
samples as well as multidirectional ‘star burst’ dispersal of 
extant haplotypes. The biogeographic modelling suggest 
three vicariance events before 5 Mya and multiple dispersal 
events in more recent times. Thus, the expectations formu- 
lated above (cf. goals 1A—C; 2) are supported. 

The observed low genetic diversity over Amazonia 
markedly differs from many other widespread Amazoni- 


and that turned out to be complexes of species with rela- 
tively small allopatric distributions (e.g. PADIAL & DE LA 
RIVA 2009, ANGULO & ICOCHEA 2010, FUNK et al. 2012, 
FOUQUET et al. 2014, 2021b, GEHARA et al. 2014, PELOSO 
et al. 2014, Rojas et al. 2018, RÉJauD et al. 2020). Notable 
in this context is the ASAP result that revealed the exist- 
ence of two lineages that could represent distinct taxa in 
our A. trivittata sampling. In this case, a small fraction of 
samples comprising the Madeira River group (from 5 lo- 
calities with 3 unique haplotypes, Fig. 2) would specifically 
differ from all other A. trivittata samples. However, this is 
not supported by the relatively low p-distances (< 0.023; 
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Figure 3. Ancestral area reconstruction of the three-striped poison frog, Ameerega trivittata, inferred in BioGeoBEARS under the 
best-fit DIVALIKE model, with: (a) the most likely biogeographic scenario, with pie charts at nodes displaying the probability of the 
ancestral estimated areas in colors and less likely estimations in white. Colors correspond to the BR as indicated on the left margin. 
Vicariance events are V1 to V3, while all other splits represent dispersal events. Group names as in text. (b) Map of Amazonia with the 
seven BRs and biogeographical stochastic mapping (BSM) with anagenetic dispersal events; numbers indicate mean range-expansion 
dispersal events and are shown with a mean of > 0.6 only. Note that only basal nodes received robust support (*), see text. 
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Supplementary Table S2) that rather suggest conspecificy 
(cf. VENCEs et al. 2005b). 

Intraspecific gene genealogies indicated a high haplo- 
type diversity (Hd = 0.966), and several identical haplo- 
types were recorded from different localities (Fig. 2b). 
These data, in combination with the relatively high ts/tv 
ratio (11.55), suggest a recent evolutionary history and in- 
terchange of populations. This is also supported by the ‘star 
burst’ pattern of the parsimony network for samples from 
western and central Amazonia (BR1, BR4, BRs), implying 
this region to be a dispersal center for our study species. 
In line with this, haplotypes associated with localities in 
the periphery of Amazonia (BR2, BR6, BR7) were geneti- 
cally fairly distant. These findings propose that A. trivittata 
invaded into these localities. In conclusion, the intraspe- 
cific gene genealogies well match the hypothesis of a recent 
within-species evolution. 

Not only the haplotype network, but also the molecu- 
lar clock and the biogeographic modelling support that the 
many younger splits in the A. trivittata phylogeny (< ca. 
4.5 Mya) are likely the result of relatively recent forward 
and backward dispersal events among BRs. Only prior to 
that, in the Late Miocene, largely between 7 and 5 Mya, vi- 
cariance was apparently a major driver for within-species 
segregation. Stochastic mapping revealed that these ear- 
ly splits have taken place in western to central Amazonia 
(BRi, BR4, BRs5, BR6). We find support for a phylogeo- 
graphic scenario that (i) with the formation of the young 
Amazon River and the retention of enormous water mass- 
es in the Late Miocene (HoorN et al. 2010, SHEPHARD et al. 
2010), A. trivittata underwent vicariance events. (ii) When 
later, around 3 Mya, water bodies burst into the Amazon 
River and diminished, A. trivittata was able to disperse be- 
tween BRs and to expand its geographic range through in- 
vasion into the periphery of the Amazon basin. 


Translation into hypothetical 
phylogeographic scenarios 


Our results are consistent with geological and landscape 
diversification processes of the Amazon basin of the last 
10 million years. During this period, central Amazonia 
changed from the lacustrine Acre System gradually into 
the modern Amazon watershed with an expansion of terra 
firme forests (HooRN et al. 2010). Enabled by the uplift of 
the Andes, the drainage of these lacustrine accumulations 
and influxes, which formed the early stages of the Ama- 
zon River, was able to breach the Purus Arch around 10.1- 
8.3 Mya (SHEPHARD et al. 2010, GORINI et al. 2014, ALBERT 
et al. 2018). The formation of the Amazon basin obvious- 
ly was a major geodispersal event with predictable conse- 
quences for biodiversity, despite that rivers were not the 
only drivers of species diversification in Amazonia (CRA- 
CRAFT et al. 2020). As figured out by Hoorn et al. (2010), 
these changes in the landscape played a decisive role in 
the evolutionary and biogeographic history of Amazonian 
plant and animal biodiversity. 


The three-striped poison frog emerged in this time 
(8.44 Mya; 95% HPD: 10.6-6.04 Mya; Fig. S3), which is not 
only supported by our time estimations but also by those of 
other authors (SANTos et al. 2009, GUILLORY et al. 2020). 
According to SANTOS et al. (2009), A. trivittata is of Am- 
azonian origin in the vicinity of the remnant of the for- 
mer Pebas System. We propose that like other contempo- 
rary biota (REJAuD et al. 2020, FouquUET et al. 20214), north 
and south of the remnant of the Acre System, it colonized 
continuous terra firme forest in central Amazonia via the 
opening of the Purus Arch (Fig. 4a, b). However, we can 
only assume this here because of vicariance events - a re- 
sult of our biogeographic modelling - that took place later. 

Since the Late Miocene (11.6-5.3 Mya), the Andes in the 
West more and more attained their present shape. Only 
since then, the formation of the Amazon basin as present- 
ly known began. The establishment of the transcontinental 
Amazon River and the arrangement of its tributaries took 
until about 4.5 Mya (GORINI et al. 2014, ALBERT et al. 2018) 
or even longer (SHEPHARD et al. 2010). Rivers that initi- 
ated through the continuous Andean uplift and today are 
known as major tributaries of the Amazon, for longer peri- 
ods were enormous water accumulations, because the Pu- 
rus Arch in central Amazonia had not entirely eroded. Our 
findings of early vicariance scenarios (Fig. 3a) suggest that 
at least three of these arising water catchments represented 
temporary riverine barriers to A. trivittata. 

One vicariance event (V1) that took place around 
6.87 Mya coincides well with river captures of the Upper 
Madeira basin (Fig. 4c) (Albert et al. 2018, Cooxe et al. 
2012, TAGLIACOLLO et al. 2015). Another vicariance event 
(V2) noted in A. trivittata occurred at around 5.75 Mya. 
It can well be linked to the emergence of the Negro Riv- 
er (Fig. 4d). This water catchment resulted from the turn 
of the Branco River, which originally drained north into 
the Caribbean Sea via the Essequibo River (ALBERT et al. 
2018). It then had become an enormous water catchment 
for about 1-2 Mya north of the Amazon River (ALMEI- 
DA-FILHO & MIRANDA 2007, RIBAS et al. 2012). The Ne- 
gro River had acted as a barrier for various taxa, including 
primates or even birds (pD’Horta et al. 2013, BOUBLI et al. 
2015). While the water retention of the young Negro River 
might have represented a riverine barrier to A. trivittata, 
we find no signal that the lower Amazon River was im- 
permeable to this species. At least, LOUGHEED et al. (1999) 
showed that for extant populations of Allobates femoralis, 
another Amazonian terra firme poison frog, the Jurua Riv- 
er is not a barrier, which perhaps is about the same size as 
was the lower Amazon River at about 6 Mya. 

The three putative isolates of A. trivittata remained with 
the complete dissolvement of the Acre System. Afterwards, 
with the forebulge of the Fitzcarrald Arch at about 5 Mya, 
geomorphology and drainage systems in the South of the 
upper Amazon basin massively re-shaped, including a shift 
of the upper Amazon (Solimões) River to the North (EspurT 
et al. 2010). The Jurua and Purus Rivers might then have 
turned into barriers to A. trivittata, corresponding to a pos- 
sible third vicariance (V3) event at around 5 Mya (Fig. 4e). 
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Guiana 
Shield 


Brazilian 
Shield 


PR = Pebas Remnant PA = Purus Arch ES = Esequibo River B = Branco River 
AS = Acre System TA = Tapajós Arch AM = Amazonas River T = Tapajos River 
VA = Vaupes Arch GA = Gurupá Arch M = Upper Madeira River J = Jurua River 
IA = Iquitos Arch FA = Fitzcarrald Arch N = Negro River P = Purus River 


Figure 4. Hypothetical biogeographic history of the three-striped poison frog, characterized by three vicariance events through river 
barriers, indicated by yellow, red, green and orange colors. Main mountains (tan), ancient arches (stipples) and water bodies (blue) 
are shown, with arrows indicating flow direction. Underlying is a map of extant northern South America. Main rivers are adopted 
from illustrations in ALBERT et al. (2018) and Hoorn et al. (2010). 
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For the last 3.5 Mya, the landscapes of the Amazon basin 
have not undergone drastic changes anymore, only rivers 
slimed down at some extent (HoorN et al. 2010) so that 
apparently A. trivittata was able to cross them and to ex- 
pand its geographic range further (e.g. into the Guianas), 
exhibiting its extant pan- Amazonian distribution (Fig. 4f). 


Dispersal advantages as a key explanation 


When trying to understand the evolutionary biogeography 
of Amazonian taxa, species’ intrinsic characteristics, includ- 
ing dispersal ability and niche breadth are pivotal (e.g. P1- 
RANI et al. 2019, FoUQUET et al. 2021b). Our data infer that 
A. trivittata once it has expanded its range into central Ama- 
zonia (similar to other anurans; e.g. SANTOS et al. 2009, RÉ- 
JAUD et al. 2020), riverine barriers had a lesser effect on this 
species than on other anurans, so that time of segregation 
was too short to evolve irreversible reproductive isolation 
(lineages). Afore mentioned species-specific characteris- 
tics may have played a key role here, and we suggest that 
the three-striped poison frog exhibits life history traits that 
make it a more successful disperser than other anurans. 

First, as mentioned above, A. trivittata is one of the larg- 
est poison frog species. Its maximum adult size is about one 
third larger than that of Allobates femoralis sensu lato (StL- 
VERSTONE 1976, LOTTERS et al. 2007), and thus, A. trivittata 
is larger than many other Amazonian anurans (see OLIVEI- 
RA et al. 2017). In anurans, larger body size favors disper- 
sal (WOLLENBERG et al. 2011, SEARCY et al. 2018). Moreover, 
larger amphibian species commonly have a higher body 
mass and can better cope with unfavorable conditions com- 
pared to smaller species (DUELLMAN & TRUEB 1986). This 
may upsurge survival rates in accidental long distance dis- 
persal over rivers known in anurans (cf. MARIN DA FONTE et 
al. 2019). Such an event, although probably rare, was report- 
ed in the Amazonian poison frog Ameerega hahneli (Mon- 
TANARIN et al. 2017). It is smaller than A. trivittata, so that 
long distance dispersal success, undergoing potential fluvia- 
tile drift, should even be higher in the larger A. trivittata. 

Dispersal success in A. trivittata is additionally stimu- 
lated by its diurnal and terrestrial life style (LOTTERs et al. 
2007, KAHN et al. 2016). Most Amazonian anurans are noc- 
turnal so that in this region mainly poison frogs cover this 
general niche; they are all smaller than A. trivittata (Mon- 
TANARIN et al. 2017). An exception is Ameerega bassleri, 
co-occurring with A. trivittata along the lower Andean 
versant in Peru (ie. A. bassleri does not emerge into low- 
lands). It is about the same size and clearly outcompetes 
A. trivittata (TwoMEy et al. 2008). In other sites, where 
various poison frog species occur in syntopy, A. trivittata is 
the most common one (SCHLUTER 2005; authors’ unpubl. 
data), suggesting that probably it better asserts itself along 
limiting niche axes. 

Other aspects increasing the dispersal ability of A. tri- 
vittata can be linked to particular behavioral traits of pa- 
rental care. Poison frogs actively promote offspring disper- 
sal by carrying their larvae to water bodies over relatively 


large distances (PaSuKONIS et al. 2019). We propose that 
the length of this transport way is positively correlated with 
adult size and that A. trivittata performs active offspring 
dispersal over larger distances than smaller poison frogs. 
At least, PaSuKoNIs et al. (2019) demonstrated that this 
species carries its tadpoles for up to 300 m away from its 
home range (comparative data for other taxa are sparse). 
Moreover, A. trivittata is one of the poison frog species 
with the highest number of eggs, probably also a result of 
its large size. It is known to lay up to about 50 eggs in one 
clutch (inferred from one male transporting 46 tadpoles 
and one female having 51 ovarian eggs; AICHINGER 1991). 
Known maximum clutch size in other Ameerega species 
is considerably less, mostly < 20 (AICHINGER 1991, LOT- 
TERS et al. 2007). Anuran clutch size is positively correlat- 
ed with dispersal ability (TRAKIMAs et al. 2016), and it has 
been shown to represent an important aspect in a disper- 
sal study on West African amphibians (PENNER & RODEL 
2019). Compared to other anurans, for instance explosive 
breeders, clutch size in A. trivittata is ‘tiny. However, in 
combination with the species’ brood care behavior includ- 
ing active offspring dispersal, dispersal via reproduction 
may be quite effective in it. 


Future research 


Poison frogs are an important group with regard to under- 
standing Amazonian biogeography (SANTos et al. 2009). 
Our results suggest that one member of this group, A. tri- 
vittata, represents an extraordinary case, because it encom- 
passes a truly pan-Amazonian distribution. Despite that 
we have an incomparable spatial sampling and provide an 
explanatory phylogeographic scenarios, at the same time 
the accuracy of our study is limited by the use of only one 
mtDNA marker and with a varying number of samples per 
locality. Moreover, some questions remain unanswered, es- 
pecially those concerning spatial processes in more recent 
times (< 3.5 Mya). The use of nuclear genomic data might 
provide a clearer picture and provide results beyond ours. 
However, the availability of materials is challenging and 
gathering nuDNA data for all the terminals included in our 
study was out of reach, so that our study is best be seen asa 
frame for a future approach (cf. Réyaup et al. 2020). Future 
research might also aim to better emphasize the dispersal 
advantages that we suggest A. trivittata has compared to 
other poison frogs. 


Conclusion 


The evolutionary biogeography of pan-Amazonian anuran 
taxa is little understood. We, for the first time, show that 
recent (Pliocene/Pleistocene) multidirectional dispersal 
plays a role, using the three-striped poison frog, A. trivitta- 
ta, as a case study. Due to its intrinsic morphological and 
life history traits it is a more successful disperser than other 
anuran species, so rivers are more permeable. 
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